O-terminated interface for thickness-insensitive transport properties of aluminum oxide Josephson junctions

Alumina Josephson junction has demonstrated a tremendous potential to realize superconducting qubits. Further progress towards scalable superconducting qubits urgently needs to be guided by novel analysis mechanisms or methods to reduce the thickness sensitivity of the junction critical current to the tunnel barrier. Here, it is first revealed that the termination mode of AlOx interface plays a crucial role in the uniformity of critical current, and we demonstrate that the O-terminated interface has the lowest resistance sensitivity to thickness. More impressively, we developed atomically structured three-dimensional models and calculated their transport properties using a combination of quantum ballistic transport theory with first-principles DFT and NEGF to examine the effects of the Al2O3 termination mode and thickness variations. This work clarifies that O-terminated interface can effectively improve the resistance uniformity of Josephson junction, offering useful guidance for increasing the yield of fixed-frequency multi-qubit quantum chips which require tight control on qubit frequency.

www.nature.com/scientificreports/ at the macroscopic level facilitates the adjustment of the process parameters to improve its homogeneity, but the absence of theoretical mechanisms at this level is not conducive to the improvement of the fabrication process. Another approach is computational modeling, which is low cost and high-throughput over the former experimental method. Simmons model was used to describe the conductance of Al/AlO x /Al junction, which did not require excessive computational resources, but the presence of empirical parameters inevitably led to errors and inaccuracies in the calculation results 16 . Jung et al. 17 calculated parameters from first-principles to fit experimental data. Koberidze et al. 18 investigated how six different Al/Al 2 O 3 interfaces affected transport properties through the potential barrier at the microscopic level by constructing an atomic model of a heterojunction. Koberidze et al. 19 showed that alumina surface irregularities and interfacial and internal relaxation affected the width and height of the barrier layer, thus altering the uniformity of the oxide layer and having a serious impact on the device performance. The transport properties of alumina tunnel junctions were simulated with non-equilibrium Green's functions (NEGF) by Cyster et al. 20 , which indicated an exponential dependence of junction resistance on oxide density. DuBois et al. 21 presented a methodology for constructing atomic-scale computational models of Josephson junctions using a combination of molecular mechanics, empirical and ab initio methods. Diešková et al. 22 performed a computational study of possible local geometric structures of interfaces at the ab-initio DFT/ GGA level of approximation to complement recent experimental data on ultra-thin AlO x -based interfaces. Even small atomic rearrangements on the AlO x -based interfaces played a significant rule in transport properties of Josephson junctions. Despite numerous efforts on the better control of the junction homogeneity both from experimental and computational perspectives, an adequate alumina interface is still absent. In comparison with the materials-centric history of classical digital computing, quantum systems have yet to discover their silicon and silicon oxide-the perfect interface between silicon and silicon oxide is the main reason for the silicon era.
Here, we reveal for the first time that the termination mode of the Al 2 O 3 interface plays a crucial role in the transport properties of the junction, and demonstrate that the conductance of the O-terminated mode has the lowest sensitivity to thickness. We developed atomically structured three-dimensional α-Al 2 O 3 models and calculated their transport properties using a combination of quantum ballistic transport theory with first-principles DFT and NEGF to examine the effects of the Al 2 O 3 termination mode and thickness variations. Furthermore, the termination mode had a more significant effect on transport properties than the thickness of similar oxide layers and the O-terminated interface can effectively improve the resistance uniformity of Josephson junction. Here, we emphasize the impact of interfacial termination mode on the transport properties of Al/AlO x /Al Josephson junctions using three-dimensional full-device models, providing valuable insights into solving the frequency deviations problem.

Results
Designing atomic-scale models of Josephson junctions with different termination modes. The atomic-scale models of Josephson junctions are based on the first-principle methods, as implemented in the Nanodcal software 23 . Previous studies have shown that matched Al(111)/α-Al 2 O 3 (0001) interfaces are the most stable structures 24,25 . Therefore, we used Al(111)/α-Al 2 O 3 (0001)/Al(111) models shown in Fig. 1, where the xand y-directions were periodic and electrical transport was along the z-direction. The structures had left and right electrodes and a central scattering region. The electrodes extended to z = ±∞ , where bias voltages were applied and electric currents were collected. The periodic structure had a semi-infinite length. Buffer layers in the boundary areas of the scattering region shielded the electrodes from scattering effects 23 , and the buffer layer material was consistent with that of the electrodes. Three different termination modes were included for the Al 2 O 3 (0001) surface 28,31,32 : O-terminated mode via one layer of O atoms, Al-terminated mode by one layer of Al atoms, and 2Al-terminated mode via two layers of Al atoms.
We also investigated the effects of various termination modes on junction transport properties for various Al 2 O 3 thicknesses. We fixed the left side of Al 2 O 3 as O-, Al-, and 2Al-terminated modes, and increased tunnel barrier thicknesses on the right side layer by layer with O, Al, 2Al-terminated modes respectively. 12 models with different thicknesses were obtained with the same termination mode on the left side of Al 2 O 3 , and 36 models with different thicknesses were calculated in total.

Zero-bias conductance for different termination models.
To qualitatively characterize the conductance differences between the device models with three different Al 2 O 3 termination modes, we first calculated the conductance without applying any bias voltage, i.e., zero-bias voltage. The zero-bias conductance for models with different interfacial termination modes is shown in Fig. 2a. Figure 2a shows that the conductance of the three models differ significantly as the oxide layer is extremely thin, which verifies that the junction conductance is very sensitive to interfacial effects as suggested by Müller et al 33 . The conductance of O-terminated model is an order of magnitude lower than the other two models, and the 2Al-terminated model has the largest conductance. These values are in good agreement with the transmission coefficients near the Fermi level in the equilibrium state, as shown in Fig. 2b, where the transmission is the greatest for the 2Al-terminated model and lowest for the O-terminated model. This may be because the Al-and 2Al-terminated models form metallic channels, enabling electrons to pass through. In contrast, the O-terminated junction model has the smallest transmission coefficient at the Fermi energy (0 eV), resulting in the smallest conductance. www.nature.com/scientificreports/  www.nature.com/scientificreports/ using Landauer formalism, as shown in Fig. 3a. The currents of all three models exhibits linearly increasing ohmic behavior when the voltage is close to 0 V, and an exponential increase at high voltages. This is consistent with experimental measurements 34 and previous calculations for thin-film tunneling 35 . Transmission spectra for the O-, Al-, and 2Al-terminated models at different voltages are shown in Fig. 3b-d, where the purple dashed lines at the left and right ends correspond to the bias window. The energy levels of Al 2 O 3 and Al are different at different biases. When the bias window is small, the transmission spectrum is flat and the current increase almost linearly. The transmission coefficients increase significantly when the bias window is large, and new peaks appear because the number of channels that participate in electron conduction is not the same under different biases and the transmission of each channel is not necessarily the same, which indicates that the current increases non-linearly with increasing biases, in agreement with the I−V curve in Fig. 3a. The transmission of 2Al-terminated model is relatively large, which corresponds to the large current obtained from the Landauer formalism.

Forward and reverse I−V curves for different termination models.
To clarify whether the asymmetry of the internal atomic structure has an important effect on the transport properties of the Josephson junctions, we calculated the forward and reverse currents of the models. Voltages in opposite polarities were applied to each of the three device models to determine differences between forward and reverse currents. The forward and reverse I−V curves for the O-, Al-, and 2Al-terminated models are shown in Fig. 4. Small differences in the forward and reverse currents are qualitatively consistent with previous measurements 17 . In the transmission spectra shown in Fig. 3b-d, the chemical potential of the whole system is different when the applied voltage is opposite in polarity. This corresponds to the integration of the spectra over a certain range on the left and right sides of the Fermi-energy level. The spectra are not perfectly symmetrical on both sides of that level because of small differences in the interface structures at the two ends of Al 2 O 3 . Therefore, the current obtained by integrating the spectrum over a certain range through Eq. 8 differs when various voltage polarities are applied. However, these small differences have no effect on the qualitative analysis of the junctions. O-terminated model has the smallest difference, while Al-terminated model has the largest difference, indicating that different termination modes affect the responses of the different models to various voltage polarities.   Fig. 5a-d. We calculated the zero-bias conductance of the models having the same termination on the right side, while having the O-, Al-and 2Al-terminated models on the left side respectively, as shown in Fig. 5e-g. Figure 5e-g shows that the zero-bias conductance decays exponentially with the thicknesses of the potential barrier layer increasing, which is in good agreement with the theory and previous reports 16 . The conductance with 2Al on the right side is larger because of metal channels and thus increased electron tunneling probabilities. For comparison, the zero-bias conductance for different thicknesses for all the 9 models are plotted in Fig. 5h. The best conductance of the device models is achieved when both ends of the oxide layers are structured with 2Al termination. The conductance is worst when both ends were O termination, and the conductance of all the models decreases exponentially because the probabilities of carrier tunneling decrease sharply. The termination mode plays a dominant role in determining the conductance of different junction models with similar thicknesses. For example, when the thickness of 2Al-terminated oxide layer is 22 Å, the conductance is larger than that of the O-terminated model having an oxide layer thickness of 20 Å. This may be because the oxide layer terminated with 2Al has metallic properties, while the oxide is an insulator in general. The termination mode of the alumina/Al interfaces has a significant effect on the transport properties, which suggests improvement possibilities in the fabrication process of Josephson junctions.

Relative change rates of conductance for models with different thicknesses. Small changes in
the oxide layer thicknesses could significantly affect the junction conductance 37 . To evaluate this sensitivity, we calculated the relative conductance change rates against the barrier thickness for different alumina termination modes at zero bias. We defined the relative conductance change rate as (G n − G n+1 )/d Å, where G n is the conductance of a junction with a Al 2 O 3 barrier of n layers, and d is fixed at 22 Å, at which the conductance is almost invariant with the increasing thickness. As shown in Fig. 6, for 9 different oxide termination modes, the change rates for each mode decreases with the increasing thickness. Noticeably, the O-O model shows a lowest change rate, while the 2Al model has the highest one. This indicates that the O-O model is more insensitive to thickness, as compared to other models.

Discussion
The structural and nanochemical properties of thin AlO x layers are decisive for the transport properties of Josephson junctions, which are crucial for qubit frequencies and the scaling of superconducting qubits. In this paper, we characterized three different Al 2 O 3 interface termination modes, and calculated their zero-bias conductance and the I-V curves accordingly. In addition, the effects of termination modes on junction transport properties for different thicknesses were also considered. The results indicate that the interfacial termination mode of Al 2 O 3 in the Al/Al 2 O 3 /Al system greatly affects the electronic transmission and the zero-bias conductance. The 2Al-terminated models have the highest conductance, while O-terminated models have the lowest conductance. The conductance is exponentially dependent on the oxide thickness regardless of the termination mode of the system. In a comparison between the thickness of Al 2 O 3 and the termination mode, the termination mode has a greater effect on the transport properties. In addition, there is a significant difference in the sensitivity of the  www.nature.com/scientificreports/ where asterisks denote surface species. The process of alumina fabrication process using ALD is shown in Fig. 7. Firstly, an atomic layer terminated by hydroxyl groups is formed on the surface of the substrate. Then trimethylaluminum (TMA) is introduced, which reacts with the surface hydroxyl groups (OH − ) to deposit a layer of Al(CH 3 ) 2 on the surface. Then, inert gas is introduced and water is hydroxylated on the surface reactant again, thus forming a complete cycle. Repeated this process can obtain a continuous and uniform high-quality Al 2 O 3 film. Based on the chemical reaction formula, it is theoretically possible to fabricate Al 2 O 3 films in the O-terminated mode by ALD, so that Josephson junctions with lower thickness sensitivity can be obtained, thus improving the resistance uniformity of Josephson junctions. However, the current experimentally fabricated Josephson junctions using ALD techniques do not have ideally O-terminated oxide layers due to inevitably formed oxide layers in between the Al films forming process and ALD in addition to other fabrication limits. Although precise control of the alumina termination modes cannot be achieved using current fabrication techniques, the schematic diagram of ALD process for fabricating oxide layers proves that it is theoretically feasible to control the termination (2) AlCH * 3 + H 2 O → AlOH * + CH 4 Figure 6. Relative change rates of conductance with increasing thicknesses for 9 termination mode models. www.nature.com/scientificreports/ mode to be O-terminated which is favorable for precise control of qubit frequencies. We expect that in the near future, fabrication process might be improved towards controlling the termination modes of the oxide layers.

Methods
The relation between qubit frequencies and transport properties of Josephson junctions. Qubit frequency follows ω 01 = ( 8E J E C − E C )/ , where E J is many times greater than E C and E J = I C 2e . The critical current I C is set by the tunnel barrier of area ~ 200 × 200 nm and 1-3 nm thickness and is thus challenging to fabricate with precision better than a few percent. According to the Ambegaokar-Baratoff relation 39 I C = π� 2eR n (where is the superconducting gap energy), normal state resistance R n relates to I C and is readily measurable to precision better than 0.1%.
The NEGF-DFT method. To calculate the transmission coefficient T(E) that an electron at energy E scatters from the left electrode (l) to the right electrode (r), we first calculated the Green's function: where G R and G A denote the retarded and advanced Green's functions of the model respectively, S is the overlap matrix due to orbital non-orthogonality (if the orbits are orthogonal, S is the unit matrix), H is the Hamiltonian matrix of the system, R and A denote the delayed and overrunning self-energy respectively. G R can be calculated from a matrix transformation using G A = (G R ) † . With the Green's functions, the transmission coefficient was calculated with: where R l (E)( A l (E)) is the retarded (advanced) self-energy of electrode l. The Brillouin zone sampling k-point was set to 20 × 20 × 1 and 100 × 100 × 1 for these calculations respectively. When the results of the two data sets were compared, the accuracy differed by only 1%, while the calculation speed was improved by a factor of twelve when the k-point was set to 20 × 20 × 1.
The conductance is given by: where G 0 = 2e 2 h is the conductance quanta, h is Planck constant, e is the electron charge, f l (E) − f r (E) is the difference between the Fermi distributions of the left and right electrodes, V r − V l is the difference between the bias applied to the left and right electrodes, and T(E) is the transmission coefficient at a certain energy E. When the bias approached zero, the limit was also the conductance at equilibrium.
The conductance is related to the transmission of the carriers, where the tunneling probability of the carriers through the one-dimensional square potential barrier is given by: where D 0 is a constant that is close to 1, a represents the barrier width (which may be different from the physical thickness of the oxide layer), and U 0 is the barrier height. Hence, the transmission coefficient decreases sharply if the barrier width or height is increased, while the barrier height is determined by the thickness and nature of the insulating layer. However, the one-dimensional square potential barrier does not completely describe our model because there are effects of various interface terminations. It is usually described by a trapezoidal potential barrier 18 , but the effects of barrier width on the conductivity are qualitatively consistent.
The current is given by the Landauer formalism: which indicates that the current is independent of how the voltage is applied, but only related to the difference noted above. The current could be considered a projection of the transmission spectrum in a certain energy range.
Quantum ballistic transport theory. Our system is a mesoscopic system between macroscopic and microscopic, which can be obtained from experimental fabrications. The length of the middle Al 2 O 3 part of the system is only a few nanometers, which is much smaller than the average free range of electron motion, when the electron transport process can be considered as quantum ballistic transport 40 . In the actual process, when the Al 2 O 3 length is short enough, the system will still be affected by temperature, material defects and other factors, which may lead to diffusive transport processes, but the ballistic transport process still plays a dominant role. Therefore, our calculations on the basis of this theory for various models can achieve a qualitative analysis of the electrical transport properties of Josephson junctions.
(3) www.nature.com/scientificreports/ Protocol. We used Nanodcal 23 first-principle quantum-transport software based on NEGF-DFT, without empirical parameters, to calculate the electrical properties of the models and to predict current-voltage characteristics and electron transmission probabilities for various junction models. During the self-consistent cycle, we used the Perdew-Burke-Ernzerhof functional in the generalized gradient approximation 41 . The convergence criteria for both the Hamiltonian and density matrices were set to 10 −5 eV, and we used the linear combination of atomic orbitals (LCAO) basis set to expand the Kohn-Sham (KS) wave function, the plane wave truncation energy was set to 100 Hartree, the maximum number of steps chosen for self-consistency and the self-consistent mixing ratio were set to 200 and 0.05 respectively, and the k-point was set to 5 × 5 × 100 and 5 × 5 × 1 respectively in self-consistent calculation of electrode part and central part.
Geometric optimization was performed with a projector-augmented wave method based on DFT [42][43][44][45] . The Perdew-Burke-Ernzerhof functional in the generalized gradient approximation was used for the exchangecorrelation interactions between electrons 41,46 . The cut-off energy of the plane-wave expansion was 400 eV, and the Brillouin zone was sampled with a Gamma-centered scheme, using 4 × 3 × 5 and 4 × 4 × 1 k-point sampling to optimize the Al and Al 2 O 3 cells. The total-energy convergence criterion was at 10 −5 eV per atom, and the maximum Hellmann-Feynman force deviation was less than 0.01 eV/Å.
Using optimized materials to build the crystal model, we varied interface contact distances and ensured that those of Al and Al 2 O 3 at the left and right ends were the same, and calculated the variations in single-point energies with the contact distance to obtain the lowest energy point. The distance corresponding to the lowest energy was the optimum distance. The optimum interface distances for O-terminated, Al-terminated, and 2Alterminated models were 2.9 Å, 3.1 Å, and 3 Å, respectively.
Several atom layers were optimized in the model where Al and Al 2 O 3 were in contact. The Brillouin zone was sampled with a 4 × 3 × 1 mesh of the Gamma-centered k-point, and the remaining parameters were the same as those in the initial material optimization. Finally, models with O, Al, and 2Al terminations shown in Fig. 1 were obtained.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.